HPC chain Stage 7 correctness: 2nd-pass Percolator, file_stems-aware reconciliation hash, Stage 6 hard-error gates - #37

The cal_match diagnostic dump previously used `{:.10}` "everywhere so we
don't hit banker's vs round-half-up rounding differences between Rust
and C#". In practice the formatters still disagree at the 10th decimal
on f64 values that land on a rounding boundary: Rust's `{:.10}` uses
round-half-to-even while .NET Framework 4.7.2's `F10` uses
round-half-away-from-zero, producing 1-in-the-last-digit (1e-10) print
diffs even though the underlying f64s are bit-equal.

Cross-impl strict comparison on Stellar 3-file 466K calibration matches:
under `:.10` the apex_rt column showed 15,739 rows "differing" by exactly
~1e-10. Under `:.17` and parsed back to f64, every matched row's apex_rt
is bit-equal cross-impl. The "drift" was 100% printf rounding artifact.

17 fractional digits is enough to round-trip any f64 uniquely, so the
diagnostic dump now exposes only real f64-level divergence and not
formatter noise. The remaining cal_match diffs (correlation, libcosine,
xcorr, snr) are real f64-level drift from the f32-throughout xcorr
preprocess on the Rust side; that is a separate concern.

Same rationale as the cal_match :.10 -> :.17 change in commit 1089407:
Rust's `{:.10}` (round-half-to-even) and .NET Framework 4.7.2's `F10`
(round-half-away-from-zero) disagree on f64 values that land on a 10th-
decimal rounding boundary by 1 in the last digit, producing a fake 1e-10
"drift" between cross-impl dumps even when the underlying f64s are
bit-equal.

Cross-impl strict comparison on Stellar 3-file 466K calibration matches:
under `:.10` the LDA discriminant/q-value showed 17 rows differing by
1e-10. Under `:.17` and parsed back to f64, max diff drops to 5.06e-14
(under 500 f64 ULP at value ~0.95) and most LDA scores are bit-equal
modulo 1 ULP.

The remaining LDA drift comes from each impl training LDA on a
different per-entry calibration-match accumulator (Stage 3 cross-impl
divergence in which apex spectrum is picked per entry) — separate
concern.

…arity

Rewrites the body of SpectralScorer::xcorr to mirror C# OspreySharp's
SpectralScorer.XcorrAtScan bit-for-bit: bin in f64 with
`(intensity as f64).sqrt()` (widen-before-sqrt, matching
`Math.Sqrt((double)float)`), apply windowing normalization and sliding-
window subtraction in f64, sum unique fragment bins in f64, scale by
`0.005_f64`. Removes the old allocating wrappers (`apply_windowing_-
normalization` and `apply_sliding_window`); the `_into` variants stay
for the HRAM per-window cache path.

Calibration apex_rt, correlation, and libcosine columns of cal_match
are now bit-equal cross-impl on Stellar Single (was 5.55e-16 / 4.97e-14
/ 5.55e-16 noise, all at f64 epsilon).

Cross-impl xcorr column drift is now 3.876e-6 (was 5.24e-10): not a
regression. The previous baseline was both impls on f32, where both
sides did pure-f32 windowing/sliding-window so the f32 cascade errors
matched bit-for-bit. With Rust flipped to f64, C# calibration -- which
still uses `XcorrFromPreprocessed(float[])` against a pure-f32 cache
in PerFileScoringTask.cs:2275 -- diverges by exactly f32 magnitude.
The coordinated fix requires C# calibration to switch to
`XcorrAtScan` (existing f64 path) for the apex xcorr; tracked in TODO.

HRAM main-search hot path (preprocess_spectrum_for_xcorr_into and the
per-window f32 cache it feeds) is untouched: memory budget preserved.

See ai/todos/active/TODO-20260516_ospreysharp_wsl_parity.md

Implements the "f64 scratch, f32 storage" architecture user-requested for
Option 3 of the cross-impl calibration parity work. Both the calibration
inline path (SpectralScorer::xcorr, surgical patch in 690194a) and the
HRAM main-search cache build path now run the windowing / sliding-window
cascade in f64 on per-thread pooled scratch, narrowing to f32 only at
the final per-spectrum cache write.

Cross-impl results on Stellar Single (joined by entry_id+charge+scan):

* cal_match: apex_rt bit-equal, xcorr drift 3.876e-6 -> 4.429e-8
  (~100x tighter). correlation, libcosine at f64 epsilon as before.
* Stage 4 .scores.parquet: peak_apex 100% bit-equal, apex_rt
  462,375/462,802 bit-equal, rt_deviation max 2.32e-11 (LOESS cascade
  now sees bit-equal inputs), xcorr max 5.41e-7 at the f32 single-cast
  architectural floor (was f32 cascade).

Changes:

* xcorr_pool.rs: XcorrScratch.binned/windowed/prefix Vec<f32> -> Vec<f64>.
  Per-thread scratch memory grows from ~1.2 MB to ~2.4 MB per worker on
  HRAM (~19 MB total at 16 threads). Per-spectrum cache (Vec<Vec<f32>>)
  unchanged at ~400 KB per spectrum.
* lib.rs apply_windowing_normalization_into: f64 in, f64 out. Mirrors
  C# ApplyWindowingNormalizationD bit-for-bit.
* lib.rs apply_sliding_window_into: f64 spectrum + f64 prefix scratch
  in, f32 result out. Single deterministic cast at final store.
* lib.rs preprocess_spectrum_for_xcorr_into: bin with
  (intensity as f64).sqrt() (widen-before-sqrt, matches
  Math.Sqrt((double)float)). Output signature unchanged (&mut [f32])
  so the HRAM per-window cache callers in pipeline.rs need no updates.
* scorer.xcorr() inline body (from 690194a) remains the canonical
  reference for the f64 cascade; it does not yet share helpers with
  the new f64 _into path (cosmetic dedup deferred).

See ai/todos/active/TODO-20260516_ospreysharp_wsl_parity.md

Closes the residual ~1e-7 cross-impl xcorr drift after Option 3 by
matching C#'s accumulator type and code path:

* xcorr_sparse: f32 accumulator -> f64 accumulator (widen each cached
  value to f64 on read, sum in f64, scale by 0.005_f64). Mirrors C#
  XcorrFromPreprocessed(float[]) which does `double xcorrRaw += preprocessed[bin]`
  (implicit float->double promote on read).
* scorer.xcorr(): the inline f64 preprocessing body is replaced by
  delegation to xcorr_at_scan, which goes through the f32 cache
  (preprocess_spectrum_for_xcorr returns the same f64-internal-narrowed
  Vec<f32> the HRAM main-search cache uses). Now bit-equal with C#
  calibration's XcorrFromPreprocessed(windowPreprocessedF32, entry)
  call site at PerFileScoringTask.cs:2275 because both impls do the
  same operations on the same f32 cache values.

Cross-impl results on Stellar Single (joined by entry_id+charge+scan):

* cal_match xcorr: 4.43e-8 -> 5.11e-15 (f64 epsilon, was f32 cast floor)
* cal_match LDA scores: 3.29e-5 -> 4.95e-14 (f64 epsilon, cascade from
  bit-equal LDA inputs)
* Stage 4 .scores.parquet xcorr: 100% bit-equal (462801/462801)
* Stage 4 .scores.parquet sg_weighted_xcorr: 100% bit-equal
* Stage 4 .scores.parquet peak_apex: 100% bit-equal
* Stage 4 .scores.parquet apex_rt: 462375 bit-equal + f64 epsilon
* rt_deviation: max 2.32e-11 (LOESS cascade, bulk <1e-12)

Only cross-impl drift remaining in cal_match is snr (5.24e-10), which
is the LDA in-place mutation hypothesis - separate root cause, next.

When the second-pass calibration refinement is accepted (R² check at
pipeline.rs:1105), Rust now also:

* Overwrites the LOESS_INPUT diagnostic dump with pass 2's points.
  Previously dump_loess_input was only called at line 1000 (before
  the refinement block) so the file reflected pass 1's 6,398 points
  even when pass 2's 7,361 points were what the LOESS fit actually
  used. C# overwrites the dump unconditionally on pass 2; this
  brings Rust into parity. Stellar Single now passes the LOESS_INPUT
  boundary bit-equal (6400 -> 7361 rows on the dump, matching the
  actual calibration data).
* Updates num_confident_peptides metadata in calibration.json to the
  refined count. Was reporting pass 1's value (6,398) even when pass
  2's 7,361 was the count actually used for the fit.

Both changes are bookkeeping/observability fixes: the rt_calibration
model parameters themselves were already bit-equal cross-impl at f64
epsilon (~1e-13) because both impls used the same 7,361 pass-2 points
for the LOESS fit. Only the diagnostic dump and the metadata count
were stuck on pass 1.

Remaining cross-impl divergence in calibration.json is ms1_calibration
{count, mean, median, sd, adjusted_tolerance} — Rust 18 errors vs C#
193. Separate root cause (ms1_error collection path).

The MS1 isotope envelope extraction at calibration time (batch.rs:2708
in run_coelution_calibration_scoring) takes a tolerance_ppm argument
and uses it for FindPeakPpm-style matching of the M+0 peak. Rust was
passing config.precursor_tolerance.tolerance regardless of unit. For
unit-resolution data (Stellar config: tolerance=1.0, unit=Mz), this
treated 1.0 Da as 1.0 ppm and produced a ~0.5 mDa window at 500 m/z -
effectively zero. As a result envelope.has_m0() returned false on
~99.8% of matches and only 18 of 7361 calibration-passing entries
contributed to ms1_calibration on Stellar Single, vs 193 on C#
(10x undercount).

C# OspreySharp's PerFileScoringTask.ScoreCalibrationEntry handles the
same case via `unit == Ppm ? tolerance : 10.0`. This adds a Rust
helper `ms1_envelope_tolerance_ppm` with identical behavior and wires
it into both pass-1 and pass-2 run_coelution_calibration_scoring
callsites in pipeline.rs.

Cross-impl effect on Stellar Single calibration.json:
ms1_calibration.{count, mean, median, sd, adjusted_tolerance} drift
disappears from the diff (was the dominant CAL_JSON divergence with
diff 1.75e+2 at ms1_calibration.count). Remaining cal_json drift is
rt_calibration.model_params.abs_residuals[i] sort-order swaps on a
handful of indices, where the underlying values are bit-equal but
sorted slightly differently cross-impl.

Found via OSPREY_DIAG_MS1 instrumentation showing tol_ppm=1.0 and
peaks_in_m0_window=0 on all sampled calibration entries.

Stable sort_by on x alone preserves input order for ties. When inputs
arrive in different orders cross-impl (e.g. discriminant-score-sorted
with 1 ULP differences), the duplicate-x positions end up at different
indices and the LOESS abs_residuals[i] diverge cross-impl. Sorting by
(x, y) makes the order deterministic and removes the input-order
dependency. C# OspreySharp LoessRegression.cs has the matching fix
(LINQ ThenBy on y).

Two small alignments bring Stage 4 .scores.parquet and Stage 5/6 dumps to
byte-identical cross-impl on Stellar net8.0 in the per-side OspreySharp
parity harness:

* compute_cosine_at_scan: switch to single-pass dot+norms with one final
  divide (matches the cosine_angle helper in the same file and the C#
  port's ComputeCosineAtScan), bringing sg_weighted_cosine to bit-equal
  cross-impl. The previous per-element divide-then-sum form is mathema-
  tically equivalent but differs in the last 1-2 ULP.

* write_scores_parquet_with_metadata: iterate entries in canonical
  (entry_id, charge, scan_number) order before filling the column
  builders. Physical row order in the parquet was previously whatever
  upstream deduplication left it in, and Rust's entry_id ascending
  pattern did not match the C# port's target-decoy-paired pattern. The
  Stage 5 standardizer and the SVM working-set selection sum / iterate
  in physical row order, so per-side cross-impl runs hit a sum-order
  cascade in the standardizer even though every column was logically
  bit-equal. Sorting before write isolates Stage 5+ from parquet-writer
  row-order differences.

Both changes are pure no-op on a single side reading itself back; only
the cross-impl behaviour changes. cargo fmt / clippy / test all pass.

The 2nd-pass `.fdr_scores.bin` was previously written only when
reconciliation was enabled across multiple files. In single-file mode
that gated the sidecar off, while the OspreySharp port writes it
unconditionally — so per-side cross-impl Stage 7 (`--join-at-pass=2`)
hit a load-cached-vs-retrain asymmetry: C# loaded the cached 2nd-pass
scores and Rust re-trained the 2nd-pass SVM from scratch, producing
different SVM weights and therefore different downstream q-values.

With the gate removed, both sides write a 2nd-pass sidecar containing
the post-compaction FDR scores (which in single-file mode equal the
1st-pass scores, since no rescore happens). The cross-impl sidecars
are then bit-identical (verified on Stellar Single, net8.0) and the
load path is symmetric.

This closes the Stage 6 → Stage 7 boundary as a source of cross-impl
drift; any remaining Stage 7 divergence is in the protein-FDR /
parsimony code itself, not in upstream inputs. Single-file workflows
gain ~5-9 MB of extra disk per file for the new sidecar; the file is
ignored by all current consumers except the resume path that wanted
it anyway.

cargo fmt / clippy / test all pass.

Cross-impl bisection of the C# OspreySharp port revealed that its
--join-at-pass=2 path was filtering second-pass detected_peptides
through stale 1st-pass q-values rather than reloading the 2nd-pass
FDR sidecar onto post-compaction stubs (Rust does this at the
pipeline.rs:4480-4494 reload block). The discrepancy was 19 peptides
on Stellar Single — bordering peptides that pass 1st-pass FDR at <=1%
but not 2nd-pass — and produced a 1-protein delta in the Stage 7
picked-protein output.

The dump is gated on OSPREY_DUMP_DETECTED_PEPTIDES=1 and writes
rust_stage7_detected_peptides.txt next to the working directory.
Zero overhead when unset. Matches the cs_stage7_detected_peptides.txt
dump the C# port writes under the same env var, so a sorted-line
diff localizes any remaining detected_peptides drift to the specific
modified_sequence values.

The 1st-pass input is already entry_id-sorted via deduplicate_pairs
(pipeline.rs:6123), but the post-rescore pool that feeds 2nd-pass
Percolator can have gap-fill entries appended after the sorted
pre-existing rows. Re-sorting at the top of run_percolator_fdr
guarantees identical iteration order across Rust and OspreySharp;
without it, gap-fill ordering can drift the cross-impl 2nd-pass
SVM working-set selection on multi-file datasets even when feature
columns are bit-equal.

On Stellar Single this is a no-op (no reconciliation → no gap-fills).
On multi-file the sort shifts the 2nd-pass dump count slightly
(Stellar 3-file: 5372 → 5360 proteins in the dump). The remaining
cross-impl drift at Stage 7 for multi-file is in code outside the
input order and is still under investigation.

cargo fmt / clippy / test all pass.

The streaming Percolator path runs a best-per-precursor dedup before
peptide-group subsampling (pipeline.rs:5512-5557) precisely so the SVM
trains on one observation per precursor instead of N observations per
N-file experiment. The direct path was originally written without that
dedup, which is statistically incorrect: on multi-file inputs sized
between max_train_size and max_train_size * 2 (e.g. Stellar 3-file at
393k entries) the SVM trained on N-times-redundant precursor pairs as
if they were independent samples. The C# port (OspreySharp.FDR.Percolator
direct path) had already inferred the correct dedup-then-subsample shape;
mirroring it here brings the two implementations into algorithmic
agreement and removes the multi-file artefact from the direct path's
training set.

Dedup key is `features[0]` (fragment_coelution_sum, the first PIN
feature) which matches the streaming path's `coelution_sum` field
value-for-value. Single-file is unaffected (each base_id has one
observation; dedup is a no-op). cargo fmt / clippy / test all pass.

…agate experiment q-values across base_id

Two coupled fixes that take 3-file Stellar from Rust 5,366 vs C# 6,541
protein groups (stage 7 FAIL) to bit-equal `.2nd-pass.fdr_scores.bin`
sidecars and stage 7 PASS.

1. **`run_percolator_fdr` sort key**: extend the per-file canonical sort
   from `(entry_id, charge, scan_number)` to
   `(entry_id, charge, scan_number, parquet_index)`. On 3-file Stellar the
   post-reconciliation pool contains 94 per-file groups of 2 entries that
   share all three of the original keys -- a gap-fill rescore landing on
   the same scan as an original row with a different `rt_deviation`. With
   only three keys, Rust's stable `sort_by` left those ties in input
   order, but .NET `List<T>.Sort` is unstable and swapped them, drifting
   the 2nd-pass standardizer mean by 1 ULP on `rt_deviation` and
   cascading through every downstream SVM weight. `parquet_index` is
   intrinsic to the byte-equal cross-impl parquet layout, so adding it
   as a final tie-break makes the total order identical on both sides
   regardless of underlying sort stability.

2. **`rescore_per_file_loop` gap-fill remap**: after rewriting the
   reconciled parquet, every `FdrEntry.parquet_index` whose row was
   moved by the writer's canonical sort must be remapped to its post-
   sort position. Previously upstream rows held a stale pre-sort index;
   when the next Percolator pass loaded features by `parquet_index`,
   the lookup silently fetched a different entry's feature row. The new
   code computes the canonical permutation explicitly (matching the
   writer's sort key), inverts it into pre->post row indices, then
   walks `fdr_entries` and remaps every stub -- both upstream rows
   (`old_pq_idx -> pre_to_post[old_pq_idx]`) and gap-fill stubs (via
   the new `gap_vec_idx_for_pre_sort_row` mapping built during append).
   Comment on `write_scores_parquet_with_metadata` documents that its
   internal sort is now stable so external pre-sort + tie-breakers
   survive the writer.

3. **`compute_experiment_precursor_qvalues` propagation**: the direct
   Percolator path assigned the winning q-value only to the single
   `compete_all` winner per `base_id`, leaving every non-winning
   per-file observation at `q=1.0`. The streaming path
   (`pipeline.rs::run_percolator_fdr_streaming`) already propagated via
   its `base_id_exp_prec_q` map; the OspreySharp port matched the
   streaming semantics. The asymmetry silently broke downstream stages
   that gate on `experiment_precursor_qvalue` (Stage 6 consensus
   selection / calibration refit, Stage 7 protein FDR) on multi-file
   inputs sized below the streaming threshold (Stellar 3-file at 393K
   entries). Direct path now builds a `base_id -> q` map from winners
   and propagates to all observations.

Verified with
`Test-Regression.ps1 -Dataset Stellar -Files All -StartStage stage6
-StopAfterStage stage7 -Tag perside_3file_v4`: stage6 PASS on all four
compare dumps (multicharge / consensus / reconciliation / rescored),
stage7 PASS on protein FDR. `diff_fdr_bin.py` shows 0 score / q-value
diffs cross-impl on every per-file `.2nd-pass.fdr_scores.bin`.

A separate one-shot diagnostic dump `dump_stage5_perc_input` (gated by
`OSPREY_DUMP_PERC_INPUT=1`) is added to localize future standardizer
divergence; writes `rust_stage5_perc_input.tsv` with per-entry raw
feature vectors sorted by `(entry_id, native_position)`.